
d=load('6.28.all');
place=d(:,1)
lati0=d(:,2)*1e-5+30.74
longi0=d(:,3)*1e-5+103.92
latitude=d(:,4)*1e-5+30.74
longitude=d(:,5)*1e-5+103.92

radiusEquatorial = 6378.1370*1000;
radiusPolar = 6356.7523*1000;

%args=[latitude longitude zeros(length(latitude),1)+latitude(1) zeros(length(longitude),1)+longitude(1)]
args=[lati0	longi0 latitude longitude]

thetaAll=arrayfun(@theta,args(:,1),args(:,2),args(:,3),args(:,4))
avg=mean(thetaAll)
result=arrayfun(@distance,args(:,1),args(:,2),args(:,3),args(:,4),zeros(length(args),1)+avg)


%lati=latitude(2:37)-latitude(1:36)
%longti=longitude(2:37)-longitude(1:36)

%args=[lati longti zeros(36,1) zeros(36,1)]
%result=arrayfun(@distance,args(:,1),args(:,2),args(:,3),args(:,4))
%result=abs(result)

plot(place,result,place,place*2.5-1.25);
result=result-place*2.5+1.25;
x=0:0.1:10;
result=abs(result);
m=mean(result);
s=std(result);
d=power(s,2);
f=(erf((x-m)/(s*sqrt(2)))+1)/2;
p=exp(power(x-m,2)/(-2*d))/(s*sqrt(2*pi));
figure
plot(x,f)
figure
plot(x,p)
